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Abstract 

Multiple  scale  methods,  which  are  based  on  discrete  and  continuous  reproducing 
kernels,  wavelets,  and  integral  window  transforms  are  developed  .  In  this  development,  a 
microscope  is  constructed  with  a  flexible  space-time  localized  window  function  which 
translates  and  dilates  in  space  and  time  to  cover  the  entire  domain  of  interest.  This 
microscope  can  magnify,  examine,  and  record  the  image  of  the  various  scales  and  frequencies 
of  the  response  locally  within  the  support  of  the  window  function.  The  degree  of  magnification 
will  depend  on  the  power  of  the  microscope,  a  flexible  space-scale  and  time-frequency  window 
function.  This  complete  characterization  of  the  unknown  response  is  performed  through  the 
integral  window  transform.  This  localization  process  can  be  achieved  by  dilating  the  flexible 
multiple-scale  window  function.  The  zoom  in  and  zoom  out  capability  of  the  window  function 
is  especially  useful  in  examining  complex  flow  phenomena,  such  as  flow  induced  vibration, 
dynamic  stability  of  flow-structure  interaction,  turbulence  structures,  and  high  frequency 
structural  dynamics  response. 
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1.  Introduction 

In  a  compressible  flow-structure  system,  the  unstable  response  often  arises  from  a 
coupling  between  phenomena  associated  with  substantially  different  frequencies.  For 
example,  vortex  formation  in  a  flow  about  a  blade  or  an  airfoil  is  initiated  by  relatively  high 
frequency  modes  of  the  structure.  The  excitation  forces  generated  by  the  vortex  and  the 
instability  itself  generally  involve  low  modes  of  structural  response.  The  complete  time 
integration  of  the  equations  for  the  compressible  fluid  and  structure  combined  can  be 
prohibitively  time-consuming  in  even  two  dimensions  and  is  beyond  the  capability  of  even  the 
largest  supercomputers  in  three  dimensions.  Thus  it  can  be  seen  that  methods  which  can 
effectively  treat  problems  with  large  ranges  in  scale  are  needed  in  many  types  of  analysis 

which  arise  in  structural  dynamics  problems. 

The  development  of  computational  methods  for  low-frequency  response  3D  structural 
systems  are  now  reasonably  well  established.  On  the  other  hand,  little  success  has  been 
made  in  the  structural  area  characterized  by  multiple  scales,  where  response  is  often 
dominated  by  the  middle  part  of  the  spectrum.  The  author  believes  that  the  multiple-scale 
reproducing  kernel  particle  methods  coupled  with  wavelets  proposed  here  have  great  promise 
in  dealing  effectively  with  the  difficult  structural  response  problems  in  which  the  medium 
frequencies  are  important,  such  as  problems  involving  impact,  dynamic  instability, 
compressible  flow-structure  interaction,  and  other  local  phenomena. 

Because  the  frequency  shift  enables  the  time  step  to  depend  only  on  the  size  of  the 
frequency  band  and  not  on  the  frequency  extent  of  the  load,  these  methods  will  speed  up 
analysis  of  medium  frequency  response  immensely.  This  is  particularly  attractive  for 
structural  response  exhibiting  multiple  time  scales.  This  new  development  will  enable 
engineers  not  only  to  bring  more  detail  into  their  structural  system  models,  but  will  also 
enhance  the  computer  simulations  of  many  classes  of  multi-scale  structural  dynamics 
analysis.  It  will  improve  the  accuracy,  efficiency,  and  reliability  of  dynamic  analysis. 

In  the  next  section,  a  review  of  the  coupled  compressible  flow-structure  interaction  is 
presented.  In  section  3,  the  weak  form  of  the  equations  of  motion  is  given  and  some  currently 
available  solution  methods  are  discussed.  In  section  4a,  the  discrete  orthogonal  reproducing 
kernel  interpolation  functions  are  reviewed.  In  section  4b,  the  multi-resolution  wavelets 
analysis  is  given.  In  section  5,  we  present  our  proposed  approach  to  study  complex 
structures,  the  multiple  scale  reproducing  kernel  particle  methods.  Sample  examples  are 
given  in  section  6. 

2.  Review  of  Coupled  Compressible  Flow-Structure  Interaction 
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The  ability  to  solve  general  classes  of  fluid-structure  interaction  problems  involving 
finite  deformations  and  stability  depends  upon  the  ability  to  solve  the  corresponding 
uncoupled  fluid  and  structural  problems,  and  also  the  ability  to  interface  fluid  and  structural 
subdomains.  During  the  last  two  decades,  considerable  progress  has  been  made  in  the 
solutions  of  free  and  moving  boundary  problems  which  involve  large  fluid  deformations. 
Among  these  methods  are  the  marker-and-cell  (MAC)  methods  (Hirt  [1975,  1983],  Nichols 
and  Hirt  [1971],  Harlow  and  Welch  [1965]),  the  volume  of  fluid  (VOF)  methods  (Amsden 
and  Harlow  [1970],  Harlow  et  al.  [1976]),  moving  mesh  techniques  (Subbiah  et  al  [1989]), 
Eulerian  and  arbitrary  Lagrangian-Eulerian  (ALE)  methods  (Liu  et  al.  [1988,  1991],  Huerta 
and  Liu  [1988],  Belytschko  [1983],  Belytschko  and  Kennedy  [1978]),  the  smoothed  particle 
hydrodynamic  (SPH)  methods  and  the  free  Lagrange  methods  (Monaghan  and  Gingold 
[1983],  Gingold  and  Monaghan  [1982],  Burton  and  Harrison  [1991]).  However,  the  coupling 
of  these  methods  with  deformable  structures  is  not  well  understood  and  often  causes 
difficulties. 

In  the  MAC  method,  a  fixed  or  Eulerian  mesh  is  used  for  the  fluid  calculation  and  a 
Lagrangian  set  of  marker  particles  is  used  to  trace  the  moving  free  surface.  Marker  particles 
can  be  spread  over  all  fluid  occupied  regions  with  each  particle  specified  to  move  with  the  fluid 
velocity  at  its  location.  A  free  surface  is  defined  as  lying  at  the  “boundary”  between  regions 
with  and  without  marker  particles.  More  specifically,  a  mesh  is  said  to  contain  a  free  surface 
if  the  mesh  contains  cells  with  markers  and  has  at  least  a  neighboring  cell  with  no  markers. 
The  MAC  method  offers  the  distinct  advantage  of  eliminating  all  logic  problems  associated 
with  intersecting  surfaces.  This  method  is  also  readily  extendable  to  three  dimensional 
computations.  However,  the  MAC  method  requires  significant  increase  in  running  time  and 
storage. 

In  the  VOF  approach,  a  function  F  is  defined  whose  value  is  unity  at  any  point  fully 
occupied  by  fluid  and  zero  otherwise.  Cells  with  F  values  between  zero  and  one  must  then 
contain  a  free  surface.  This  fractional  volume  of  fluid  (VOF)  method  provides  the  same 
coarse  interface  information  available  in  the  marker  particle  method.  In  order  to  trace  the  free 
surface,  a  flow  analysis  network  approach  is  developed.  It  uses  the  flux  calculation  at  the 
boundary  of  each  control  volume  to  check  the  fraction  of  fill  instead  of  tracking  marker 
particles.  One  limitation  of  this  control  volume  technique  is  that  the  time  step  must  be 
controlled  to  ensure  that  free  surface  only  passes  one  control  volume  in  a  time  step. 

In  a  moving  mesh  approach,  a  numerical  grid  generation  scheme,  which  facilitates 
solutions  over  arbitrarily  shaped  boundaries,  has  to  be  developed.  This  approach  numerically 
maps  the  irregular  shape  of  the  flow  field  to  a  more  regular  shape  in  a  computational  domain 
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where  the  governing  equations  are  solved.  The  computation  time  of  this  method  is  relatively 
high  because  of  the  large  number  of  time  steps  and  the  mesh  generation  at  each  time  step. 

There  arc  two  approaches  to  ALE  methods.  One  approach  updates  the  solution 
variables  in  a  single  time  step  while  the  other  performs  a  Lagrangian  step  follows  by  an 
Eulerian  or  remap  (or  advection)  step.  The  latter  strategy  is  often  referred  to  as  an  operator 
split  method.  Similar  to  the  moving  mesh  approach,  a  new  mesh  is  required  for  the  ALE 
calculations.  However,  this  new  mesh  is  required  only  if  the  Lagrangian  mesh  is  too 
distorted.  The  major  advantage  of  the  ALE  methods  is  the  moving  boundaries  can  be 
computed  with  high  accuracy  of  the  Lagrangian  method,  and  the  mesh  can  conserve  its 
regularity  in  avoiding  element  entanglement  Although  a  combination  of  multi-material  ALE 
methods  have  also  been  developed,  the  application  of  this  technique  to  multiple  free  surfaces, 
multiple  materials,  and  fluid-structure  interaction  remains  to  be  explored. 

The  SPH  (smooth  particle  hydrodynamics)  and  the  free  Lagrange  methods  are  based 
on  algorithms  which  are  truly  grid-free  or  mesh-free.  ConsequenUy,  these  techniques  do  not 
have  the  mesh  entanglement  problems  and  at  the  same  time,  maintain  the  high  accuracy  of 
the  Lagrangian  calculations.  Similar  to  a  finite  difference/element  Lagrangian  calculation,  the 
continuous  fluid  is  approximated  by  a  set  of  particles  or  fluid  elements  of  equal  mass.  Unlike 
the  traditional  Lagrangian  finite  element/difference  calculations,  the  inter-particle  forces 
among  these  smooth  fluid  particles  are  derived  from  the  pressure  and  these  interacting 
pressure  forces  are  governed  by  an  interpolating  function.  Hence,  the  nodal  or  element  forces 
(which  are  usually  constructed  from  a  finite  element/difference  mesh)  are  no  longer  needed; 
and  consequently,  these  free  Lagrange  methods  require  no  mesh  and  the  mesh  distortion 
problems  are  completely  eliminated.  The  energy  equation  is  similarly  defined  for  each  smooth 
fluid  particle.  Since  the  fluid  is  modeled  by  equal  mass  particles,  the  density  for  the  fluid  is 
constructed  by  defining  a  weighting  function  in  which  the  density  of  the  equal  mass  particles 
is  proportional  to  the  number  of  particles  per  unit  volume.  With  these  approximations,  the 
motion  of  the  fluid  is  governed  by  the  movement  of  this  set  of  smooth  fluid  particles,  and  the 
movement  of  the  particles  are  governed  by  the  particle  interacting  forces. 

Although  SPH  methods  work  well  if  there  is  no  boundary  (since  the  boundary  terms 
are  tossed  out  in  the  formulation,  Libersky  and  Petschek  [1990]),  and  when  the  number  of 
unknowns  (nodes)  is  large;  SPH  methods  are  not  as  accurate  as  the  regular  finite  element 
methods,  Johnson,  Peterson  and  Stryrk,  [1993].  From  our  study  of  SPH  interpolation  function 
via  a  simple  one  dimensional  (ID)  Galerkin  formulation,  we  found  that  there  is  an  additional 
deficiency  in  the  SPH  formulation.  It  is  related  to  the  boundary  correction  term  of  the 
reproducing  kernel  approximation.  We  shall  make  an  attempt  to  identify  this  deficiency  and 
present  our  view  of  improving  the  SPH  kernel  approximation. 
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After  reviewing  the  moving  least  square  interpolation  functions,  Lancaster  and 
Salkauskas  [1981],  and  the  diffuse  element  methods  (DEM),  Nayroles  et  al.  [1992]. 
Belytschko  et  al.  [1993]  pointed  out  that  an  assumption  made  by  Nayroles  et  al.,  the 
interpolation  coefficients  are  constants,  detracts  from  the  accuracy  of  the  method.  They 
developed  the  Element  Free  Galerkin  Methods  (EFGM)  and  showed  that  by  adding  more 
accurate  derivatives  and  enforcing  boundary  conditions  by  Lagrange  multipliers,  the  methods 
could  achieve  very  high  rates  of  convergence.  From  our  experience,  EFGM  are  more  accurate 
than  the  finite  element  methods,  and  hence,  the  SPH  methods  especially  for  a  small  set  of 
nodes.  One  main  drawback  of  EFGM  is  the  computational  expense,  and  we  found  that  it  is 
more  computationally  intensive  than  the  SPH  methods. 

The  objective  of  the  Reproducing  Kernel  Particle  Methods  developed  by  Liu  et  al. 
[1993],  is  along  the  same  line  of  development  as  the  SPH,  DEM,  and  EFGM  :  to  develop  an 
accurate  and  efficient  mesh  free  interpolation  functions.  A  detailed  discussion  on  smooth 
particle  methods,  the  diffuse  element  methods  and  the  element  free  Galerkin  method,  and  the 
recently  developed  reproducing  kernel  particle  methods,  is  also  given  in  the  paper. 

Since  a  continuous  reproducing  kernel  can  be  derived  for  this  method  and  it  is  also  a 
free  Lagrange  particle  method,  we  shall  label  this  development  as  Reproducing  Kernel 
Particle  Methods  (RKPM).  This  proposed  approach  is  motivated  by  the  theory  of  wavelets 
Chui  [1992]  and  Daubechies  [1992],  in  which  a  function  is  represented  by  a  combination  of 
the  dilation  and  translation  of  a  single  wavelet,  which  is  a  window  function.  In  a  wavelet 
analysis,  similar  to  the  SPH  interpolation  kernel,  the  interpolation  coefficients  are  defined  in 
terms  of  the  integral  window  transform  of  the  window  function  and  the  solution  itself.  In  this 
proposed  study,  we  shall  make  use  of  the  multiple  frequency  bands  and/or  multiple  scales 
properties  of  wavelets  analysis  (multi-resolution  analysis),  and  the  time- frequency  and/or 
space-scale  localization  properties  of  the  continuous  and  discrete  reproducing  kernel 
approximations. 

3.  Weak  Form  of  the  Equation  of  Motion  and  Solution  Methods 

Consid^  the  weak  form  of  the  equations  of  motion  which  can  be  written  as: 

K  <5uh,  u*‘>  +  B<5ub(xerg),  u*>(xErg)  -  g(x)>  =  f<5uh,  p>  (3.1) 

where  K<-,  •>,  B<*,  •>  and  f<-,  •>  are  the  usual  weak  form  operator  of  the  governing 
equations,  boundary  constraint  operator  on  fg,  and  the  force  assembly  operator,  respectively. 
5u^(x,  t),  u^(x,  t),  p(x,  t)  and  g(x,  t)  are  the  test  function,  trial  function,  body  force,  and 
prescribed  data  on  the  boundary  fg.  The  governing  equations  can  be  those  of  structural 
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dynamics,  fluid  dynamics,  coupled-fluid-structure  interaction  or  structural  acoustics,  and 
among  others.  The  classical  Galerkin  method  is  to  approximate  u^  (which  is  the 
approximation  to  u(x.  t))  by: 

A 

u^Cx,  t)  =  X  C.  <l>(x  -  x.) 

ft—  I 


and 


A 

bu^X.  t)  =  X  5Ca  <|)(X  -  Xa)  ^b) 

where  (t>(x  -*a)  can  be  the  global  finite  element  shape  functions.  Hughes  [1987],  spectral 
functions  Gottlieb  and  Orszag  [1977],  the  smooth  particle  hydrodynamic  (SPH)  interpolation 
kernel  function,  Lucy  [1977]  and  Monaghan  [1988],  multiple  scale  finite  element  functions 
Liu,  Zhang  and  Ramirez  [1991],  or  wavelet-type  bases,  Chui  [1992]  and  Daubechies  [1992], 

etc.  Substitute  Eqs.  (3.2)  into  Eq.  (3.1)  and  solve  for  Ca  for  a  =  1 . A  coefficients  wiU  result 

in  the  discrete  approximation  of  u(x),  provided  certain  continuity  requirements  are  met 

It  is  noted  that  the  above  solution  procedures  also  hold  for  the  various  approaches  such 
as  the  space-time  discontinuous  finite  elements,  Hughes  and  Hulbert  [1988],  and  Shakib  and 
Hughes  [1991],  deforming  space-time  discontinuous  finite  elements,  Tezduyar,  Behr  and  Liou 
[1992];  the  arbitrary  Lagrangian-Eulerian  (ALE)  and  Eulerian-type  finite  element  methods, 
Liu  et  al.  [1991],  and  among  others.  All  these  methods  employ  the  same  type  of 
interpolations,  Eqs.  (3.2),  except  a  new  set  of  motion,  mesh  motion,  is  introduced  through 
similar  equations  (3.2).  This  additional  motion  is  used  to  control  the  mesh  or  grid 
deformation  so  that  mesh  distortion  can  be  minimized  and  the  nodal  points  connectivity 
through  the  mesh  or  grid  description  would  not  give  negative  Jacobians,  Liu  et  al.  [1988]. 

Equations  (3.2)  give  a  good  approximation  to  u(x,  t)  when  the  measure  of  the  spacing 
among  Xa,  call  it  h,  becomes  smaller  and  smaller.  One  can  then  visualize  the  coefficient  Ca  as 
an  average  value  of  u(x,  t),  and  <|>(x  -  Xa)  is  the  weighting  function.  If  h  approaches  to  zero  Ca 
approaches  to  the  true  solution  at  Xa.  However,  in  practice,  h  does  not  approach  zero  and  if 
u(x,  t)  is  a  very  nonlinear  function,  solving  for  the  average  values,  as  all  of  these  methods  do, 
would  probably  leave  out  the  fine  details  of  the  response  u(x,  t). 

Therefore,  our  objective  is,  instead  of  solving  for  the  averaged  u(xa)  (i.e.,  Ca),  we  shall 
develop  an  alternative  form  of  Eqs.  (3.2),  which  is  based  on  discrete  reproducing  kernels  and 


7 


J 


integral  window  transforms.  In  this  approach,  if  we  include  the  time  dimension  into  x,  we  can 
think  of  ^(x  -  Xg)  as  a  flexible  space-time  window  function  located  at  Xg.  Since  u(x)  is  an 
unknown  function,  our  goal  then  is  to  interpret  the  coefficient  Cg  as  a  microscope  which 
magnifies,  examines,  and  records  the  image  of  the  response  u(x)  around  Xg.  This  can  be 
achieved  by  employing  a  space-time  localized  window  function  that  can  recover  the  various 
scales  and  frequencies  of  the  response  u(x)  locally  around  Xg.  This  leads  us  to  the  use  of 
scaling  functions  and  wavelets  which  are  discussed  in  the  next  two  sections. 

4a.  Discrete  Orthogonal  Reproducing  Kernel  Interpolation  Functions 

If  <|>(x  -  Xg)  is  chosen  to  be  the  scaling  function,  Chui  [1992],  Cg  is  identified  as  the 
integral  window  transform  of  the  unknown  response  u  and  4»(x  -  Xg)  over  the  domain.  That  is, 
the  unknown  coefficients  are  to  be  constructed  so  that: 

C,  =  C,(u,  Xa)  S  <U,  <►,>  ss  I  u(x)  (j)(x-Xa)  dx 

Jv  (4.1) 

where  <u,  4t>g>  is  the  integral  window  transform,  and  V  is  the  domain  of  interest.  The 
reconstruction  formula  Eq.  (3.2a)  becomes: 

A 

u‘'(x,  t)  =  ]£  <U,  <j)a>  <|)(x-Xa) 

-I  (4.2) 


Equation  (4.2)  is  typical  for  a  discrete  reproducing  kernel  Hilbert  space.  In  seeking  for  the 
solution  uh(x,  t)  using  Eq.  (4.2)  and  a  similar  equation  for  Suh(x),  it  is  necessary  to  define 
nodal  point  xj;  nodal  uj  m  u(xj);  and  particles  with  nodal  mass  AMj  and  nodal  density  pj. 
Consequently,  the  nodal  volume  AVj  is  determined  by  AMj/pj  =  AVj  for  J  =  1, ...,  NP,  where 
NP  is  the  total  number  of  particles  inside  V.  Hence, 

NP  NP 

AVj  *  V  ;  and  ^  AMj  =  total  mass  (4.3) 

j=i  j=i 

Using  numerical  integration  in  Eqs.  (4.2)  and  the  above  definitions,  uh(x,  t)  becomes: 

A 

u‘‘(x,  t)  =  X  <l>a>  <I»(X-Xa) 

asl 


( 
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rueB(i,) 

=  £  I  «(X,-Xa)AV,  U,(t) 

•-l'-  ' 

where  xj  are  the  quadrature  points,  and  B(xa)  is  the  support  of  <j>(x-xa).  Equation  (4.4)  is  a 
discrete  reproducing  kernel  approximation  (DRKA),  since  u**(x,  t)  is  interpolated  via  (t>(x-xa) 
through  the  integral  window  transform  of  u(x,  t).  An  interesting  interpretation  of  Eq.  (4.4)  is 
as  follows.  There  are  a  total  of  A  sampling  windows  spread  over  the  domain  V.  At  each 
sampling  location  Xa,  a  microscope,  which  is  constructed  from  the  integral  window  transform 
<u,  <(>a>,  examines  the  motion  of  particles  AMj  carrying  nodal  values  uj  passing  through  the 
support  B(xa).  If  the  motion  is  very  nonlinear,  the  number  of  ”free  Lagrange  particles  AMj 
passing  through  each  sampling  window  <|)(x-Xa)  is  different  from  time  to  time.  Therefore,  Eq. 
(4.4)  maintains  the  attributes  of  the  free  Lagrange  methods,  and  with  the  appropriate  choice 
of  <|)(x-xa),  gives  more  accurate  results. 

Upon  examining  Eq.  (4.4),  since  only  a  set  of  nodal  points  or  particles  are  involved  in 
DRKA  methods,  similar  to  SPH  methods,  DRKA  methods  have  no  problems  associated  with 
mesh/grid  distortion  or  negative  Jacobian  of  the  elements  resulting  from  large  deformation. 
Hence,  it  is  suitable  for  large  deformation  and  high  velocity  flow  problems.  However,  unlike 
SPH  methods  in  which  <l>(x-xa)  in  Eq.  (3.2a)  is  the  interpolation  kernel  function  that  mimic  a 
Dirac  Delta  function  and  Ca  is  the  product  of  the  nodal  mass  AMa  and  the  nodal  value  of  the 
response  Ua  (-u(xa))  divided  by  the  nodal  density  pa,  (Monaghan  [1988]).  That  is 

u‘‘(x,  t)  =  £  <1>(X-X.)  Ua  S  2  <|>(X-Xa)  AVaUa 

a=l  a=l  (4.5) 

Comparing  Eqs.  (4.4)  and  (4.5),  if  only  one  point  integration  is  used  to  integrate  the  integral 
window  transform  <u,  ^a>»  which  is  a  very  bad  approximation  to  the  integral,  especially  for 
very  nonlinear  or  complex  u(x,  t),  the  SPH  and  DRKA  methods  take  a  similar  form.  However, 
if  more  than  one  point  quadratures  are  employed  to  evaluate  the  integral  window  transform 
<u,  <t>a>.  the  difference  between  the  two  methods  is  apparent  It  is  thus  clear  from  Eq.  (4.4) 
why  we  call  Ca  a  microscope  in  which  the  magnification  power  will  depend  on  the  number  of 
particles  examined  under  the  support  B(xa). 

4b.  Multi'Resolution  Wavelets  Analysis 


<|>(X-Xa) 


(4.4) 
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In  this  section,  we  wish  to  construct  and  interpret  Ca  as  a  flexible  power  microscope. 
This  microscope  can  zoom  in  (sharpen  the  window  function)  to  pick  up  the  detailed  structures 
of  the  response  tt(x);  and  zoom  out  (widen  the  window  function)  if  no  further  magnification  of 
the  response  is  necessary.  This  zoom  in  and  zoom  out  capability  is  especially  useful  in 
examining  complex  flow  phenomena,  such  as  flow  induced  vibrations;  dynamic  stability  of 
flow-structure  interaction,  turbulence  structures,  structural  acoustics,  high  frequency  structure 
dynamic  response,  strain  localization  problems,  and  other  engineering  disciplines. 

In  this  development,  since  Ca  is  defined  in  terms  of  the  unknown  response  u(x)  and  the 
multiple  scale  window  function  (|>(x-Xa),  Eq.  (4.2)  has  to  be  redefined  for  the  multi-resolution 
wavelets  analysis.  The  window  functions  <|>(x-Xa)  will  be  replaced  by  flexible  power  window 

functions  Vina(x);  whereas  the  integral  window  transform  coefficients  Ca  =  <u,  ^a>.  a  =  1 . 

A,  become  Cma,  which  are  defined  by  <u,  Vma>-  The  power  of  the  window  function  is 
controlled  by  an  index  m  =  0,  1, ...,  M.  The  multi-resolution  analysis  is  performed  simply  by 
the  dilation  of  the  window  functions.  That  is,  the  discrete  reproducing  kernel  formula  Eq. 
(4.2)  can  be  re-defined  as: 

u'‘(x)  =  2  X  Vn».>  Viii«(x)  =  X  X  CmaVmafx) 

m-O  a=l  “  m=0  t»l  (4.6a) 

where  the  dilation  of  Vina(x)  is  defined  through  the  integer  index  m,  and  AO  >  A1  >  ....  >  AM. 
The  arbitrary  scale,  ao  >  0,  is  usually  set  equal  to  2.  Hence,  the  flexible  space-time  and 
scale-frequency  window  functions  are  defined  by: 

¥ma(x)  =  ao“^  v(ao”(x-x  J)  ao  >  0  (4.6b) 

and  the  integral  window  transform  <u,  \)fiiia>  is  given  by: 

Cma  =  <u»  -  I  Viiia(x)  u(x)  dx 

jv  (4.6c) 

where  V  is  the  space  (or  space-time)  domain  of  interest 

As  can  be  seen  from  Eq.  (4.6b),  the  mother  wavelet  v(>t)  can  be  obtained  by  setting  m  = 
0  and  Xa=  0.  Because  of  our  choice  of  indexing,  m  =  0  corresponds  to  the  smallest  window 
width.  The  integral  window  transform  <u,  Voa^’  microscope,  can  pick  up  the  very  fine 

scale  and/or  high  frequencies  up  to  the  arbitrary  scale  ao.  For  m  =  M,  it  corresponds  to  the 
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largest  window  width,  and  CMa  =  <u.  ¥Ma>«  microscope,  can  pick  up  the  large  scale 
and/or  low  frequency  response. 

It  is  apparent  from  Eqs.  (4.6)  that  the  microscope  Cma  located  at  Xa  can  examine  the 
finest  scales  of  u(x)  simply  by  dilation.  The  magnification  factor  is  determined  by  the  mother 
window  function,  and  the  flexible  space-scale  and  time-frequency  localized  integral  window 
transform  <u,  Vma>- 

Using  a  numerical  quadrature  integration  scheme  to  discretize  Eq.  (4.6c)  gives: 
xjeB"(x.) 

C„a  =  <U,  <l^aul>  =  Z  ¥nu(Xi)AV,U, 

where  B™(xa)  is  the  support  of  the  m^  scale  window  function  located  at  Xa  and  AVj  is  the 
nodal  volume  evaluated  at  xj.  Substituting  Eq.  (4.7a)  into  Eq.  (4.6a),  we  obtain  our  desired 
discrete  multi-resolution  wavelet  analysis  of  u(x): 

M  An 

u**(x)  =  X  X  \  ^ 

in=0  a=l  '  J 

Similar  equation  can  also  be  written  for  5ul^(x).  Equation  (4.8)  can  be  written  in  a  more 
familiar  form  of  multi-resolution  analysis: 


uj 


(4.8) 


M 

uKx)  =  X  ■■■ 

nt^ 


(4.9a) 


It  can  be  seen  from  Eq.  (4.7a)  that  u''(x)  is  a  direct  sum  of  the  M-scale  solution.  For  each 
refinement  m,  the  mH>.scale  interpolation,  which  also  represents  the  m‘l‘-frequency/wave 
number  band  of  the  solution,  is  defined  by: 

Am 

uSi(x)  =  X  «Sa(x)  =  U^l(x)  +  «m2(*)  +  •"  +  U^(x) 

(4.9b) 


From  Eq.  (4.9b),  we  can  view  each  i4a(x)  as  a  flexible  power  microscope  examining  u(x) 
around  Xa  and  each  u&a(x)  is  interpolated  through  the  window  function  so  that: 
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(4.9c) 


Xj€B"(X.) 

U&a(x)  =  X 

J 

and  the  m***  scale-a***  wavelet  shape  function  (NmaiCx))  of  particle  J  is  given  by  (no  sum  on  m, 
a  and  J): 

Nmaj(x)  =  Vnia  (Xj)  AVj  Vma  (x)  (4.9d) 


From  Eq.  (4.9c),  the  summation  on  J  for  xj  e  B®(xa)  will  deHne  the  global  nodal  connectivity. 
However,  unlike  the  usual  finite  element,  the  sequence  of  numbering  uj  would  not  cause  any 
mesh  distortion  or  negative  Jacobian  problems  since  there  is  no  element 


If  additional  particles  are  added  to  the  domain  (this  can  easily  be  done  by  splitting  up 
particles),  the  adaptive  multi-resolution  analysis  can  be  summarized  as  follows: 

M 

u'Kx)  =  X  u^x)  multi-resolution  analysis 

in=o  (4.10a) 


M  /Am 

-ill 

in^  Va>l 


u^  (x)  summation  of  all  frequency  bands  of  interest 


(4.10b) 


M  /Am 

SI 

mpO  Vast 


X  Nmaj(x)  Uj 
J 


Uj  are  examined  by  the  m*  scale-a* 
wavelet  under  B“(xJ 


(4.10c) 


xjeB“(x.)  I  M  Am  \ 

X  { X  X  ^Vma(Xj)AVj\|rnia(x))  [uj  global  interpolation  functions 
J  ImpO  asl  I  (4.10d) 


As  can  be  seen,  there  is  no  change  in  the  multi-resolution  analysis  Eqs.  (4.10a)  and  (4.10b), 
since  the  dilation  index  m  and  the  number  of  window  functions  Wmg(x)  are  the  same.  There  is 
also  no  change  in  the  order  of  the  window  function  as  Vma(x)  is  derived  from  \Koa(x).  The 
only  change  is  in  the  summation  on  J  under  the  support  B>°(xa)  and  the  resulting 
” nodal/particle”  matrix  will  be  larger  because  of  the  additional  global  interpolation  functions  of 
those  nodes/particles.  Similar  procedures  can  also  be  developed  for  the  deletion  of  particles. 
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If  we  can  construct  discrete  multi-resolution  wavelet  functions  according  to  Eqs.  (4. 10), 
we  shall  have  a  truly  mesh  or  grid-free  adaptive  method  and  the  interpolation  is  defined 
through  a  set  of  arbitrary-spaced  nodal  points  or  particles.  Moreover,  if  v(x)  is  chosen  to  be 
Ck  smooth  functions,  i.e.,  smooth  in  its  k***  derivatives,  the  discrete  multi-resolution 
reproducing  kernel  interpolation  functions  are  also 

5.  Multiple  Scale  Reproducing  Kernel  Particle  Methods 

The  multi-resolution  wavelet  analysis  given  in  sections  4  is  based  on  an  infinite  domain 
assumption.  To  apply  wavelets  to  analyzing  complex  structures,  this  restrictive  assumption 
is  no  longer  valid.  Another  intrinsic  deficiency  of  orthogonal  wavelet  is  the  stringent 
requirement: 

I  x“'v(x)dx  =  0  m  =  0,  1, ...,  q 

I 

where  q  is  the  degree  of  polynomials  of  the  mother  wavelet  v(x).  From  Eq.  (5.1)  it  follows 
that  a  q^  order  wavelet  can  not  represent  1,  x,  x^,  ...,  x**,  parts  of  the  solution.  To  remedy 
these  two  restrictions,  we  shall  represent  a  function  u(x)  by: 


u(x)  =  u'^(x)  +  P(x)c 


(5.2) 


where  u'^(x)  is  the  part  of  the  solution  that  is  obtained  by  the  multi-resolution  wavelet 
reconstruction  as  given  in  section  4b.  P(x)=  (Pi(x),  P2(x), ....  Pn(x)}  and  c  =  {ci,  C2, ...» Cn}^ 
are  the  vectors  of  the  n  linear  independent  functions  and  unknown  coefficients,  respectively. 
A  superscript  T  denotes  the  transpose.  We  can  consider  the  P(x)c  term  as  the  residual 
representation  of  u(x)  within  a  bounded  domain.  Following  the  procedures  of  deriving  the 
reproducing  kernel  particle  interpolation  functions  (see  Liu  et  al.[l993]),  the  multiple  scale 
reproducing  kernel  interpolation  function  coupled  with  wavelets  can  be  shown  to  be; 


u(x)  =  u'*'(x)  +  I  C(x,  y,  ao)  ao‘  u(y)  dy 

-j*  C(x,  y,  ao)  ai,‘  u'*'(y)  dy 


(5.3) 
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where  the  parameter  ao  determines  the  size  of  the  scaling  function  <|)(x).  The  form  of  the 
correction  function  C(x,  y,  ao)  will  depend  on  the  choice  of  P(x).  It  is  noted  that  the  second 
term  in  Eq.(5.3)  is  the  reproducing  kernel  approximation  described  in  Appendix  A.  The  first 
term  is  the  multi-resolution  wavelet  part,  whereas,  the  third  term  connects  the  two 
reproducing  methods.  It  is  interesting  to  point  out  that  by  a  proper  choice  of  ao,  the 
contribution  of  the  coupling  term  can  be  shown  to  be  negligible.  With  this  construction,  the 
wavelet  and  the  reproducing  kernel  terms  give  the  high  and  low  frequency  (or  the  fine  and 
coarse  scales)  representations  of  the  solution  u.  It  is  also  noted  that  u^'Cx)  can  be  expressed 
by  other  continuous  or  discrete  multiple  scale  reproducing  kernels. 

To  further  examine  this  multiple  frequency/wave  number  bands  wavelet  approximation, 
we  let 


M  An 

U'^(x)  =  X  X  Vma>  ¥ma(x) 


nM)  t*l 


(5.4) 


Substituting  Eq.(5.4)  into  Eq.  (5.3)  yields: 

M  /m  f  f  \ 

u(x)  =  X  X  (  I  “(y)  ¥ma(y)  dy|[Vn,a(x)-¥ina(x)] 
n*=0  •=!  \Jy  } 

+  1  C(x,  y,  ao)  ao‘  <|>(^)  u(y)  dy 

Jv  (5.5) 


where  the  definition  of: 


<u,  ¥ma>  =  I  u(y)  Vnu(y)  dy 

Jy  (5.6a) 

has  been  used  in  Eq.(5.5).  The  approximation  of  the  wavelet  functions  \Kma(x)  through  the 
reproducing  kernel,  denoted  by  ¥ma(x)  is: 

¥ina(x)  =  I  C(x,  y,  ao)  <t>(^)¥ma(y)  dy 

Jy  (5.6b) 
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It  is  now  clear  from  Eq.  (5.6b)  that  Vm.(x)  is  simply  an  approximation  of  Vma(x)  via  the 
reproducing  kernel  reconstruction.  It  is  expected  that  ¥m.(x)  is  very  close  to  Vma(x)  for  low 
frequency/wave  number  wavelets.  Consequently,  the  contribution  from  the  low 
frequency/wave  number  wavelets  is  close  to  zero.  However,  depending  on  the  choice  of  ao,  <}>. 
and  v;  the  reproducing  kernel  might  not  be  able  to  reconstruct  the  high  frequency/wave 
number  part  of  the  solution  (second  term  in  Eq.  (5.5));  and  these  high  frequency/wave  number 
components  can  be  readily  picked  up  by  multi-resolution  wavelets  analysis. 

Presently,  we  are  working  on  the  theoretical  analysis  of  this  type  of  reproducing  kernel 
methods.  Upon  the  understanding  of  Eqs.  (5.6).  we  shall  implement  the  correct  <D.  and  V  into 
the  continuous  multiple-scale  frequency  bands  approximation  of  the  response: 

uh(x)  =  j  C(x.  y.  ao,  a*")  ab‘  u(y)  dy  low  frequency  band 

+  X  I  C(x,  y.ao,a“‘)a3"v(^)u(y)dy  high  frequency  band  (5.7) 

m=l  Jw 

Unlike  orthogonal  wavelets,  equation  (5.7)  holds  for  arbitrary  domains.  Discretization  of  Eq. 
(5.7)  gives  the  desired  multi-grid/multi-resolution  analysis  of  the  complex  dynamic  systems: 

NPO  -  _v 

uKx)  =  £  [C(x,  xj,  ao,  a”)  Axj  uj]  low  frequency  band 

j-i 

scaling  factor  coefficients 
NPl  _v 

+  £  [C(x,  xj.  ao,  a”)  Axj  uj]  a-i  \|/(^)  1st  higher  frequency  band 

j=i 

Ist-scale  wavelet  coefficients 


+  £  [C(x,  XJ,  ao.  a”)  Axj  uj]  a-M  M*  higher  frequency  band 

J=i 

Mth-scale  wavelet  coefficients 


(5.8) 
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It  is  noted  that  NPM  ^  Npl  ^  npO;  and  it  constitutes  an  unstructured  multi-grid  analysis. 

We  also  wish  to  emphasize  that  orthogonal  wavelets  are  not  necessary;  hence  there  is  a 
larger  class  of  window  functions  which  can  give  good  time  and  frequency  localization. 

6.  Numerical  Examples 

The  steady-state  advection-diffusion  problem  can  be  stated  for  the  one-dimensional  case 
as  follows: 

u,xx-  =  b(x)  in  il 

with  the  boundary  conditions 

u(xg)  =ui  on^g 

u,x(xh)  =  Uj  onTh 

where  is  the  domain  (i.e.,  0  £  x  ^  L  ),  Tg  is  the  boundary  within  essential  boundary 
conditions  and  Ft  the  boundary  with  natural  boundary  conditions.  The  whole  boundary  is 
r  =  FguFh  and  rgnrii={0}.  (  ),x  denotes  derivatives  with  respect  to  x.  u  is  the  scalar 
unknown,  a  given  constant,  and  b(x)  a  given  source  term.  The  parameter  a  is  the  advetive 
velocity  divided  by  the  diffusion  coefficient  The  Peclet  number  is  therefore 

Pe=a^  (6.2c) 

The  boundary  conditions  are  given.  This  problem  can  be  viewed  as  a  heat  transfer  problem 
with  convective  and  diffusive  heat  transfer.  The  source  term  b(x)  can  be  caused  by  a 
chemical  reaction.  Following  Hughes  et  al.  [9]  the  weak  form  of  equation  (6.1)  with  the  least 
square  term  can  be  written  as 


(6.1) 


(6.2a) 

(6.2b) 
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(6.3) 


f  w(u,xx-a  u.x-b)  dQ  +  I  (w,xx-a  w,x)  x  (u,xx-a  u,x-b)  dft  =  0 
in  in 


where  w  is  an  arbitrary  test  function  and  t  is  a  parameter.  Using  the  approximations  u**  and 
wi*  for  the  functions  u  and  w,  we  obtain  the  usual  matrix  equation. 

The  following  numerical  examples  use  a  highly  irregular  source  term  to  cause  a 
nonlinear  solution  with  two  peaks.  The  source  term  consists  of  two  terms  that  are  very 
similar  and  summed  together.  Each  part  is 

.  _  2ciC2eF^*~*”)  l)  SecKci(x-xo))^ 

■  (e2c:(x.xo)+l)2 

2c?  Sect(ci(x-xo))^  Tant(ci(x-xo)) 

{eCiX-Xo)+ffCiX-Xo)^ 

(^^ix-xo)+ccix-xo))  *  (l-Tanh(ci(x-xo))) 

(e^3(*-*o)u|.e-cj(x-xo))2 

*  (l-Tanh(ci(x-xo))) 

^eCj(x-xo)L(.erC5(x-xo))3 

( c^sc<^-xoi.e<ix.xo))  (i-TanKci(x-xo)))  Ci  Sech(ci(x-xo)p  \ 

(eca(x-xo)+erca(*-*o))2  (ec:(x-xoi+.e^^*-xo))  j 

where  the  parameter  xq  governs  the  position  of  the  peak,  ci  controls  the  sharpness  on  the 
right  side  of  the  peak  and  C2  controls  the  decay  on  the  left  side.  In  the  previous  equation  the 
second  indice  for  the  paramters  xq,  ci  and  C2  is  omitted.  Therefore  x©  =  xoi,  ci  =  cu  and 
C2  =  C2i  for  i=l,2. 

The  resulting  source  terra  is 

b(x)  =  kibi(x)  +  k2b2(x)  (6.5) 
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Figure  6.1  Source  term  used  in  the  examples 


with  the  parameters 


i 

1 

2 

XOi 

5.0 

2.5 

Cli 

50 

50 

C2i 

1.0 

3.0 

_ki _ 

0.5 

1.0 

where  the  parameter  k  determines  the  size  of  the  peak.  The  source  term  using  these 
coefficients  is  shown  in  figure  6.1. 

The  homogeneous  solution  for  the  advection-diffusion  equation  on  the  domain  0  ^  x  ^  6 
with  boundary  conditions  u(0)s0,  u(6)  =  1  becomes 

u^(x)  =  ^-'4  (6.6) 

e6«- 1 


where  the  particular  solution  with  b(x)  given  in  (6.17) 
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u^x)  =  kiU?(x>fk2U^x) 


(6.7) 


where 


p.  ,  l-TanKci(x-xo)) 


for  i  =  1,2 


(6.8) 


with  the  appropriate  constants  from  the  table  above.  The  parameter  a  is  set  to  a  1000  in  the 
following  examples. 

The  results  for  the  diffusion  equation  were  obtained  by  setting  the  parameter  a  to  zero 

and  setting  the  boundary  conditions  to  u(0)=0,  u(6)  =  0. 

The  results  for  the  Reproducing  Kernel  Method  are  obtained  by  using  a  Window  function 

W  of  the  form 


W(z)=  e-*' 


(6.9) 


(6.10) 


To  show  the  influence  of  the  parameter  j,  the  solution  of  the  differential  equation  is  calculated 
for  several  j.  Note  that  for  j=-2.2  the  solution  becomes  unstable  and  that  for  a  ;<0  the 
solution  tq)proaches  the  Hnite  element  solution. 

Note  that  the  wavelets  are  scaled  with  a  factor  Ax  to  scale  the  mother  wavelet  to  the  size  of 
the  mesh. 

The  solutions  of  the  Diffusion  Equation  and  the  Advection  Diffusion  Equation  using  several  different 
methods  are  shown  in  Figures  6.1  -  6.6  An  error  plot  for  the  Diffusion  Equation  is  shown  in  Figure  6.5.  The 
solution  for  31  nodes  is  not  representative,  the  source  term  is  underintegrated. 
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Figure  6.5a  Error  of  numerical  solutions  for 
the  Diffusion  Equation 
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Figure  6.5b  Error  in  the  derivative  of  numerical  solutions  for 
the  Diffusion  Equation 


31 

61 

121 

241 

0.19545 

0.18977 

0.07468 

0.01827 

0.19511 

0.18954 

0.07829 

0.06641 

RKPMi  =  *1.0 

0.19589 

0.19054 

0.07507 

0.01839 

BSSSBRH 

0.22884 

0.11999 

0.03746 

0.00463 

BCmiHBIII 

0.20792 

0.16593 

0.07138 

IRffRliEI 

0.19507 

0.19137 

0.07507 

Table  6. la  Error  of  numerical  solutions  for 
the  Diffusion  Equation 


num.  of  nodes 

31 

61 

121 

241 

Usual  FEM 

0.80144 

0.79926 

0.56435 

0.27424 

0.80144 

0.79910 

0.56383 

0.27401 

0.80283 

0.80029 

0.56555 

0.27422 

RKPMi=  1.0 

0.79341 

0.63022 

0.33317 

0.07116 

RSISlISlil 

0.79900 

0.73682 

0.54145 

BBWWBIW 

0.80149 

0.79854 

0.55935 

Table  6.1b  Error  in  the  derivative  of  numerical  solutions  for 
the  Diffusion  Equation 
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Figure  6.6  Solution  of  the  Advection-Diffiision  Equation 


using  Gaussian  Reproducing  Kernel  with  j=  -2.2 
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